Vortex-lattice formation and melting in a nonrotating Bose-Einstein condensate 
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Numerical simulations of the interference of a three-way segmented nonrotating Bose-Einstein condensate 
reveal the production of a honeycomb vortex lattice containing significant numbers of vortices and antivortices. 
If confined within a trap, the lattice subsequently melts, exhibiting a rich assortment of vortex-antivortex in- 
teractions. In contrast with nonlinear vortex production mechanisms previously described for Bose-Einstein 
condensates, the process here is shown to be primarily one of linear superposition, with initial vortex locations 
approximately described by a linear theory of wave packet interference. 
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I. INTRODUCTION 

The production of vortices has attracted great interest in 
the study of Bose-Einstein condensates (BECs) (see, e.g., 
IH 12 HI Si)- Typically, production has been through bulk ro- 
tation of the condensate cloud, such as with a "laser spoon" 
or by laser phase imprinting. These rotating systems form 
an Abrikosov lattice [5] of vortices with hexagonal symme- 
try. In contrast with the rotating BEC, in which the number 
of vortices is governed by the net angular momentum of the 
system, the nonrotating BEC can also give rise to vortices due 
to the reconciliation of initial random phase variations via the 
Kibble-Zurek mechanism |@] . 

Interference of two nonrotating BEC pieces with a repulsive 
nonlinearity, has also been shown to give rise to vortices [7]. 
In this system, analogous to the Young's two-pinhole inter- 
ferometer, the interference fringes — also known as dark stripe 
solitons — decay via a "snake instability" into a string of vor- 
tices This vortex formation mechanism relies on the 
nonlinearity of the BEC self interaction. 

Recently Scherer et al. [10] performed an experiment in 
which vortices were observed as a result of the interference of 
a three-way segmented BEC. An oblate spheroidal BEC was 
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FIG. 1 : (Color online) A laser-illuminated mask separates a pancake- 
shaped BEC, formed in an asymmetric trap, into three pieces. Upon 
removal of the illumination — and optionally the trap — the pieces in- 
terfere, forming a lattice of vortices and antivortices. 
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formed in an asymmetric trap partitioned into three regions by 
shining a laser light sheet on the condensate (see Fig.Q]). The 
partition walls were then removed at varying rates and vor- 
tices were sought after different elapsed times, during which 
the three pieces were allowed to expand and interfere. Both 
2D and 3D numerical simulations of this experiment have also 
been performed by Carretero-Gonzalez et al. Il ill . 

In this paper, we demonstrate with numerical simulations 
that vortices are produced by a three-segment BEC devoid of 
initial phase variations and show that this mechanism is, in 
contrast with the two-piece case, predicted by a linear theory, 
the development of which is related to previous work on the 
three-pinhole Young's interferometer IU2I1 . In contrast with 
Carretero-Gonzalez et al. Il ill , we have not sought to repli- 
cate the aforementioned experiment exactly. By instead fo- 
cusing attention on those initial states characterized by a con- 
stant phase, we aim to demonstrate that initial phase variations 
are unnecessary for vortex formation. In addition, we demon- 
strate interference and vortex production in the absence of a 
confining transverse trap, thereby reducing the requirement 
for the nonlinear processes at play in two-fragment conden- 
sate interference, and providing deeper insight into the argu- 
ment that the vortex generation mechanism for three symmet- 
rically arranged, well-separated pieces is, by contrast, primar- 
ily a linear process. 

By increasing the intensity of the light sheet, we have also 
been able to generate a lattice comprised of significant num- 
bers of vortices and antivortices. In the trapped system, the 
regular vortex-antivortex lattice subsequently melts, exhibit- 
ing a diversity of interactions by, for example, self-propelling 
vortex-antivortex dipoles (VDs) OH [U El Qj, H, 
rotating vortex tripoles and quadrupoles IU5lll7ll . These inter- 
actions between vortices and vortex clusters in the condensate 
include dipole scattering and annihilation events. The large 
vortex population produced in the trapped system makes it an 
excellent environment for the study of vortex dynamics, VDs 
and vortex-sound interactions |l|] . 

The remainder of this paper is structured as follows. We 
describe the production of vortices via the phasor approach 
in Sec. [II] In Sec. [Ill] we describe the numerical BEC model, 
focusing on the time-dependent and ground state models, the 
reduction from 3D to 2D, and a vorticity metric for the order 
parameter field. We present the simulation results in Sec. HyJ 
highlighting the effects of light-sheet intensity, different light 
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FIG. 2: Representing the sources with a phasor diagram, vortices 
lie at locations corresponding to a closed loop of phasors. (a) Three 
sources linearly superposed, (b) Two sources and a nonlinear inter- 
action term, (c) Four sources linearly superposed or three sources 
plus a nonlinear interaction term. 



sheet geometries and phase variation between the condensate 
pieces in the trapped and untrapped cases, whilst making con- 
nections to the linear theory. We provide a visualization of 
the lattice formation and melting, and in Sec. [V] describe the 
rich vortex-antivortex dynamics that arise after the lattice has 
melted. Finally, we summarize our findings in Sec. IVII 



II. PHASOR DESCRIPTION OF VORTEX PRODUCTION 

In weakly nonlinear systems, typical of experimentally re- 
alistic BEC systems, we expect to observe phenomena ap- 
proximately predicted by the linear theory. The creation of 
a vortex-antivortex lattice is one such observation, whose ge- 
ometry is described in the Appendix. We include this short 
section to explain the production of vortices by linear interfer- 
ence, from which the lattice is formed. 

In order to form vortices by a linear interference process, 
at least three source waves are required lfl9tl . This may be 
understood with reference to the phasor diagram Fig. [2f a) in 
which each phasor represents a source wave. The probability 
amplitude is zero at a vortex core, corresponding to a closed 
loop comprised in general of three or more phasors. Although 
two phasors may sum to zero, such a case corresponds to a 
domain-wall defect rather than a vortex. Vortex formation 
from linear superposition of higher numbers of sources has 
also been demonstrated B20I1 . corresponding to a loop of four 
or more phasors. 

In the following, the individual terms of a sum of complex 
scalar waves *Pi + ¥2 + • • • are represented graphically as pha- 
sors p; +P2 + - • • inFig.|2] Figure^ a) is a representation of the 
linear interference of three sources ¥\ + ¥2 + ^3 evaluated at 
the vortex location. Alternatively, Fig. |2|b) represents a non- 
linear system in which the third contribution is now associ- 
ated with the interaction term in the sum *Pj + ¥2 + f(^i, ^2)- 
The production of vortices by the snake instability might be 
represented by the latter picture. Thus, a loop of n phasors 
may correspond either to the interference of n sources in a 
linear system, or to n - 1 sources plus a nonlinear interaction 
term. For example, Fig.[2jc) may represent either four sources 
x ¥\ +^2+^3 +¥4 in the case of a linear system, or three sources 
plus a nonlinear interaction term *P] + ¥2 + ¥3 + f(*Pi ,^2,^3) 
in the case of a nonlinear system. 

Recently, the results of the interference of 2D condensates 
sliced into two and four pieces have been studied 112 ill . Exam- 
ination of these results for the case of four pieces suggest that 
sufficient contributions from at least three of the four source 



pieces have combined in the manner just described, which we 
refer to as a primarily linear process. For the system we study 
in this paper, we present further evidence for the role of the 
linear interference mechanism in Sec. [TV] 



III. NUMERICAL MODEL 

In this section we describe the derivation and parameters 
of our numerical model for the interference of a three-way 
segmented nonrotating BEC. We model an experiment similar 
to that conducted by Scherer et al. ifToll . 

We note the following features that distinguish our sim- 
ulation from the experiment and its modeling by Carretero- 
Gonzalez et al. Ill ill . First, in establishing the initial condition 
prior to time evolution we obtain the ground state with a uni- 
form global phase, so that no initial phase variation is allowed 
within or between the three pieces, in order to demonstrate 
vortex production as predicted by the linear theory. In con- 
trast, Carretero-Gonzalez et al. iflltl apply phase variations 
within their numerical scheme to obtain their initial condi- 
tion. Although we subsequently also apply an example of 
different relative phases to the three pieces following estab- 
lishment of the ground state, this is presented to demonstrate 
how the linear theory applies to this case. The second distin- 
guishing feature is the instantaneous removal of the light sheet 
at t = 0. In contrast, Scherer et al. B10I1 remove the light sheet 
at varying finite rates, but report a maximum efficiency of vor- 
tex generation for the maximum observable rate of removal. 
Thirdly, we exclusively model the system in 2D. As such, our 
simulation results are most pertinent to pancake-shaped con- 
densate clouds, in which the ratio of the axial to transverse 
trap frequencies is large in the sense outlined below. Because 
non-pancake condensates admit the existence of vortex lines 
with more complex geometries, investigation of these systems 
would have to be performed separately to assess the applica- 
bility of the results reported here. A final difference is that 
we present simulations with and without the transverse trap 
to demonstrate that lack of transverse confinement is not an 
impediment to vortex production, consistent with predictions 
from the linear theory. 

The production of vortices in pancake-shaped condensates 
may be numerically modeled in two rather than three spatial 
dimensions. In this case the evolution of the macroscopic 
wavefunction or order parameter ¥ of the BEC is consid- 
ered to be governed by the Gross-Pitaevskii equation (GPE) 
in(2+l)D, 
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where we have adopted plane-polar coordinates for space (r, 6) 
and t denotes time. The GPE takes the form of the (2+l)D 
time-dependent Schrodinger equation with an additional non- 
linear self-interaction term jNUq I*!*! 2 *P for a condensate con- 
taining N atoms, each of mass m. Here Uq = AnH 2 a/m de- 
pends on the s-wave scattering length a, which in this case is 
positive in order to permit the BEC pieces to expand and inter- 
fere. The multiplication factor y arises as a result of reduction 
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from the full 3D description to 2D in the manner described 
below. In this form of the GPE, the normalization condition is 
2n J q °° \¥(r,9)\ 2 rdr = 1 (see, e.g., j^SB). 

The axial dimension of the BEC is determined by the axial 
trap component mco 2 z 2 /2, whose angular frequency u> z pro- 
vides a characteristic length a z = [h/(mco z )] 1 ^ 2 . When a z 
is much greater than the scattering length a, the multiplica- 
tion factor y may be found by separating the 3D order pa- 
rameter into a product of independent axial and transverse 
components ^^(r, 6,z) ~ <5>o(zy¥(r, 9) and evaluating the 
expectation value of the transverse component *P(r, 9) along 
the axial direction z l23l l25ll . The axial condensate pro- 
file is assumed to be the well-known ground-state solution 
O (z) = (7Ta 2 y l/4 exp[-z 2 /(2a 2 )] to the ID time independent 
Schrodinger equation. The 2D Laplacian operator arises from 
separation of the full 3D Laplacian into its transverse and axial 
components, V 2 = V 2 + d\. The final result is a (2+l)D GPE 
with a nonlinear multiplication factor y = [mtoj (2nh)] 1 ^ 2 that 
represents the averaged value of the 3D nonlinear factor over 
all z- An alternative approach, and that presented here, is 
to directly evaluate the nonlinearity occurring at z = 0, cor- 
responding to the location of the peak density of the axial 
profile, and assign this value to y. Since the nonlinearity 
achieves a maximum in this plane, it arguably represents a 
good choice for modeling the vortex dynamics. In this case 
y = |3>o(0)| 2 = [mtoj (nti)] ' 2 . We have performed simula- 
tions with both approaches and find no qualitative difference 
with any of the effects reported here. The reader preferring the 
former approach may note that the factor of V2 by which the 
y factors differ may be absorbed as a change in the number of 
condensate atoms N. 

The BEC ground state is determined in the presence of a 
trapping potential V(r, 8) formed from the combination of a 
harmonic trap and light-sheet potential V(r, 9) = mcL> 2 r 2 /2 + 
L(r, 9), where u> is the angular frequency of the radial trap 
component. By assuming that the light sheet L(r, 9) is instan- 
taneously removed at t — 0, the light-sheet potential term ap- 
pears only in the time independent GPE used to determine the 
ground state. The harmonic trap term is included in the time 
dependent GPE for the trapped condensate simulations, but it 
is removed for the untrapped simulations. The parameter val- 
ues used in our simulations were a = 5.77 x 10~ 9 m for 87 Rb 
|H, 7Y = 2.6 x 10 5 ,w = 46.5 rad s" 1 and u, = 88.6rad s -1 . 

We define a vorticity metric for the order parameter field *P 
according to 

vorticity = |Vx j| dx dy, (2) 

where the integral represents integration over the numerical 
field of volume V, and j is the probability current j = ( V FV V I'- 
WP*)S/(2m). 

This metric is similar to the usual vorticity measure in flu- 
ids co — V x v, of the velocity field v = j/p where p = | v f| 2 
is the local probability density. By not dividing the probabil- 
ity current by p, extra weighting is conferred upon vortices 
located within a locally increased density over vortices in re- 
gions of lower density. The metric sums the modulus of the 
local measure over the whole field. 



IV. RESULTS 

In Fig.|3]we present simulation snapshots for the segmented 
BEC showing the initial state, the vortex-antivortex lattice 
that forms after the three pieces have interfered, and later 
stages that, in the trapped cases, follow the disintegration of 
the lattice. For corresponding movies, see 12611 . The intensity 
of the laser light sheet relative to the trap is represented by 
a dimensionless quantity /(>, where the trap potential contour 
plots (shown to the left of the time series) indicate the magni- 
tude relative to the trap. We have chosen to present results for 
Iq = 0.3 and Iq = 0.8, as these values allow exposition of the 
different behaviors of poorly and well separated BEC pieces, 
respectively. The walls are characterized, as shown in the trap 
profile plots, by wall height ho, expressed in units of energy, 
and full width at half maximum width wq, Iq = 0.3 and 0.8 
correspond to ho ~ ksXlQ nK and ks x 180 nK, respectively; 
both higher than the kg x 26 nK barrier produced by Scherer 
et al. iflOTI . Our wall width wo ~ 9 /urn may be challenging to 
produce experimentally. A comparison value is not given in 
Scherer et al. II X Oil - although the condensate profiles are indica- 
tive of wider walls. The upper time series show the evolution 
of the BEC beginning from the global ground state. Probabil- 
ity density and phase are also shown in Fig. [3] The expansion 
and interference occurs both within a confining harmonic trap 
[Figs. Oa) and (c)] following instantaneous removal of the 
light sheet, and in the absence of a transverse confining trap 
[Figs. [3jb) and (d)]. In the latter case, the condensate would 
expand beyond a finite simulation region. To allow for this 
in the numerical model, a damping term is added to the time 
evolution equation to absorb the outward propagating matter. 
Since these are 2D simulations, and as such are applicable to 
pancake-shaped condensates, experimental realization would 
most likely require maintenance of the axial trap component 
in both cases, not just to allow for establishment of the lattice, 
but also for any free-expansion stage prior to imaging. 

Scherer et al. [10] observed vortices consistent with pro- 
duction by the Kibble-Zurek mechanism. They reported that 
10% of nonsegmented condensates contain vortices. In our 
simulations, by starting from a ground state, no vortices may 
be produced by this mechanism. 

Previous work [12] has shown that the linear interference 
of three expanding monochromatic spherical waves generates 
a distorted honeycomb vortex-antivortex (VA) lattice. In the 
Appendix, we present a related linear theory for the case of 
three Gaussian wave packets, evolving in (2+l)D, clarifying 
the effect of source phase variation on the predicted vortex lo- 
cations. In this case, an infinite, regular honeycomb VA lattice 
is formed, with a Gaussian probability density envelope. The 
formation of vortices by a three-piece BEC may be understood 
as arising from the same mechanism, albeit now in a highly 
nonlinear system. The linear theory applies most directly to 
the untrapped system. For observing interaction dynamics, 
the trap presence must be maintained. However, if the inter- 
ference is instead performed with the transverse trap switched 
off, experimental measurement of the position of any central 
vortex with respect to the center and lattice parameter should 
allow determination of the phases of the initial BEC pieces to 
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FIG. 3: (Color online) Amplitude (upper) and associated phase plots (lower) for the ground-state trap potential shown as a contour-plot at 
left, which is a 2D harmonic trap with a superposed three-way light sheet of dimensionless amplitude Iq = 0.8. The trap profiles identifying 
wall height ho and width wo are taken along the dotted paths, (a) Three-segment BEC interference with a harmonic confining trap showing 
progression from the ground state, through lattice formation to a late-stage characterized by complex vortex-antivortex (VA) dynamics, (b) 
Honeycomb lattice formation and free expansion in the absence of a confining trap, (c-d) The initial piece phases are rotated to 0, 2n/3 and 
-2tt/3. The effect is to shift the lattice so that an antivortex coincides with the trap center. 

ory instead describes an equivalent effect, which predicts the 
presence of a central vortex as resulting from translation of 
the lattice as a whole. Rohrl et al. 112711 state that the pieces 
are "virtually degenerate", allowing them to be treated as co- 
herent pieces whose relative phases may vary randomly. For 
a lower intensity light sheet, tunneling through the light-sheet 
walls ensures that the phases of the separate pieces remain 
coupled. The interference pattern is said to become "locked". 
In our case, the lattice translation becomes locked in an equiv- 
alent sense to give the result in Figs.[3la) and (b). To simulate 
the effect of decoupling between the pieces, we apply global 
phase factors by rotating the phase of the BEC pieces, follow- 
ing establishment of the ground state but prior to time evolu- 
tion. 

In Figs. [3jc-d), the relative phases of two of the regions 
have been rotated from to 2n/3 and -2jt/3, as shown in the 
t — s phase plot. The t — 53 ms and t = 107 ms cases 
without a transverse trap show a vortex at the center. This 
may be understood by evaluating the predicted vortex loca- 
tions using the linear theory. To illustrate this, Fig.|4]provides 
linear simulation results for the equal [Figs. Hfa-b)] and ro- 
tated [Fig. |4]c)] phase cases for the model described in the 
Appendix. For these examples, the source positions r2 and 
f3, and the momentum uncertainty Ap were determined by fit- 
ting Gaussian profiles to the leading (innermost) edges of the 
probability density of the three BEC pieces at t = for the 
Iq - 0.8 case, i.e. from the first frame in Fig. |3ja). Because 
the model assumes that the BEC pieces are well described by 




FIG. 4: (Color online) Three linearly superposed wave packets ly- 
ing at the corners of an equilateral triangle [cf. Figs.[3tb) and (d)]. 
The amplitude, with analytically determined vortices (dark) and an- 
tivortices (light) overlaid [see Eqs. l lAlOt and l lAHH . and phase of 
Eq. dA2t are shown, (a) Equal-phase wave packets shown at ms, 
after interference at 50 ms (with and without the overlaid lattice), 
and (b) after 100 ms. (c) Phases 0, 2^/3 and -2n/3 cause a lattice 
translation. 



within a global phase factor. Figures |3b) and (d) show the 
formation of an extremely regular lattice, which compares fa- 
vorably with that formed by linear superposition in Fig. |4] 

If the three BEC pieces are sufficiently isolated from each 
other, they may acquire random relative phases Il27ll . Scherer 
et al. H 1 Of! describe the presence or absence of a central vortex 
according to reconciliation of these phases. The linear the- 
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FIG. 5: (Color online) Amplitude [and associated phase plots for (c)] for light sheets characterized by the trapping potentials shown as 
adjacent contour plots, (a-c) When compared with Fig.(3] the lower intensity light sheet (I = 0.3) allows vortices to propagate outwards with 
the condensate matter. In (a), the central region of the lattice is left devoid of vortices. In the trapped (b) and untrapped (c) cases — in which 
the relative phases of the initial pieces are set to 0, 2n/3 and -2n/3 — the central vortex is clearly visible, (d-e) Two alternative arrangements 
with (d) propeller- shaped pieces, and (e) pieces at three corners of a square are discussed in the main text. 



circularly-symmetric Gaussian pieces, which is not the case 
here, this approach was found to be better than fitting to the 
whole BEC-pieces. However, because of this departure, il- 
lustrated by the poor match of the lattice scale to that of the 
nonlinear results, the linear model is best applied qualitatively. 
Although the nonlinear dynamics of the BEC do not allow the 
analytically predicted vortex locations to be mapped directly 
onto the nonlinear simulation results, the linear theory nev- 
ertheless provides useful predictions of the generation of the 
central vortex and the lattice symmetry. 

If the phases are such that a vortex is created sufficiently 
close to the center of the trap, in the case where few other vor- 
tices are produced, such as when the wall heights are low [as 
in Figs. EJa-d)], the vortex may migrate to take up residence 
in the center of the trap. Figure [5jb) illustrates this behav- 
ior. Since flows associated with phase gradients cancel there, 
any resident vortex occupies a privileged position and is al- 
lowed to remain there until perturbed by a local change, such 
as might be caused by a passing vortex set in motion following 
melting of the lattice. 

The vorticity metric is plotted against time for four differ- 
ent initial conditions in Fig. |6{a). Figures |6jb-c) are time- 
series plots of the vortices in the side-on condensate slices. 
In these figures, period (i) denotes the initial formation of 
the regular vortex-antivortex lattice. In the subsequent pe- 
riod (ii), the regular lattice melts, as the condensate matter 
washes back towards the center due to confinement by the har- 
monic trap. During this period, the high population of vortices 
and antivortices promotes a high interaction rate characterized 
by a rich variety of vortex interaction dynamics. VA anni- 
hilation reduces the vortex population to a level where this 
rate slows and dynamics involving VDs begin to dominate in 
the final period (iii). The washing in and out is roughly cir- 



cularly symmetric, cyclically increasing and decreasing the 
probability density in the center. These global oscillations are 
visible in Fig. |6ja) due to the vorticity measure giving extra 
weight to vortices embedded in a locally increased conden- 
sate probability density. The predicted frequency of the ra- 
dial breathing mode of a harmonically trapped 2D conden- 
sate is twice the trap frequency a> fl28ll (In our simulations 
(i) = 46.5 rad s _1 ). The measured frequency from the lower 
three plot series, whose condensates have centroids coincid- 
ing with the trap center, is 93.6 + 0.6 rad s , in agreement 
with theory. From Figs. |6fb) and (c) we see that the vortices 
themselves are carried in and out by the bulk motion of the 
condensate matter, presumably changing the vortex interac- 
tion rate in turn. The interaction dynamics are discussed in 
more detail in Sec. [V] 

A lower wall height results in the pieces starting closer to- 
gether. In the linear theory in 11211 . correspondingly smaller 
r values generate smaller numbers of vortices. This is clearly 
illustrated by the smaller vorticity values in Fig. [6} a) and the 
reduced population of vortices in Fig. |6jc) over Fig. |6fb). In 
FigS-HJa-b), we see that the matter expanding from the center 
carries most of the vortices to the outer parts of the trap as the 
lattice is forming, where they may be lost altogether, further 
reducing the vortex population. 

The propeller shape in Fig. [3d) is included to show an alter- 
native mask, which also produces BEC pieces at 120° angles 
to each other, but which pushes the centers of the condensate 
pieces further apart. Accordingly, higher numbers of vortices 
are predicted by the BEC created from this mask. Indeed, this 
is the case, as can be seen by comparing the bottom two re- 
sults in Fig.|6ja), which plot vorticity for the same light sheet 
intensity Iq = 0.3 for the two 120° alternatives. 

In the linear theory in lfl2tl . the maximum number of vor- 
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FIG. 6: (Color online) Vorticity plots for selected BEC segmentation 
schemes with a harmonic trap. Small arrows along the time axis cor- 
respond to the times selected in Figs.[3]and[5] (a) Vorticity measure 
[Eq. l[2j] for I = 0.3 and 0.8 and for sources at three corners of an 
equilateral triangle, a square and propeller shaped schemes (see main 
text), (b) Vortices (dark) and antivortices (light) generated from the 
equilateral triangle scheme with Iq = 0.8. Regions i, ii and iii are 
described in the main text, (c) Fewer vortices are produced with a 
lower intensity I = 0.3 light sheet. 



tices is predicted for a source arrangement in which the pieces 
are initially arranged at three corners of a square. This pre- 
diction guided investigation of a new source configuration, 
shown in Fig. |3e), which may be realized by a cross-shaped 
light sheet mask with one quadrant open to the passage of 
light from the illuminating laser. Confirmation of the predic- 
tion is apparent by comparing the top two results in Fig.|6ja), 
which are for the same light-sheet intensity. In fact, the spac- 
ing of the condensate centroids for the 90° arrangement is 
smaller than for the 120° case. Thus the increase in vorticity 
due to the angular arrangement more than offsets the decrease 
in vorticity-generating capacity caused by reducing the piece 
spacing. Whilst we have shown a geometry with an increased 
capacity to generate vortices, it is possible that nonlinearity 
causes the absolute maximum to be achieved for another an- 
gle close to 90°. Evidence of nonlinear effects is visible in 
this case: the snake instability HI [H] is visible as a curvature 
of some fringes of the second frame in Fig.|3e). On close in- 
spection we see that the linear vortex generation mechanism 
dominates in this region, as disturbances from the third BEC 
piece propagate through and interfere with these fringes prior 
to the snake instability evolving to form vortices directly. The 
initial displacement of the condensate centroid from the center 
of the trap results in the merged condensate oscillating back 
and forth along the initial mirror symmetry plane, in addition 
to the radial oscillation mode. In a previous investigation the 



ability to suppress the snake instability was reported in a sys- 
tem containing a ring dark soliton (RDS) — a dark stripe soli- 
ton formed into a closed loop — that shared the circular sym- 
metry of the confining harmonic trap If29ll . Depending on the 
RDS parameters, this structure could preferentially decay by 
emission of radiation instead of the snake instability. In the 
three-segment BEC system, the intervention of linear interfer- 
ence can be seen as another means for suppressing the onset 
of the snake instability in dark stripe solitons. 

A possible explanation for Scherer et al. IToll not observing 
the production of large numbers of vortices along with a reg- 
ular lattice in experiments may be their application of light- 
sheet intensities below our lower intensity results. Blurring 
and possible nonuniformity of the light-sheet walls is related 
to the size and manufacturing process used for the mask — 
attention to the design of the mask optics and optical path may 
permit the creation of more uniform, sharper walls, leading to 
the lattice structures we see in our simulations. Also, real 
experimental initial conditions with a non-zero-temperature 
chemical potential introduce some phase randomness which 
would disrupt the lattice. Finally (as noted in 1 1 1]), if the con- 
densate is more spheroidal than pancake-like, the probability 
of observing vortices in the projected 2D image is reduced due 
to the extra freedom of possible vortex line geometries. 



V. VORTEX DYNAMICS 

In this section we focus attention on the striking interac- 
tions seen to occur between those vortices and antivortices re- 
maining in the condensate following the disintegration of the 
lattice. The vortices and antivortices are propelled and evolve 
within a seething background sea of sound waves. Due to the 
complexity of these interactions, we provide a phenomeno- 
logical description of just some of the rich dynamics, making 
observations of structures previously described by other in- 
vestigators. A fuller appreciation of the complexity will be 
gained by viewing the accompanying movies 1 26] . 

The dynamics described in this section will be exhibited 
primarily in 2D pancake-shaped condensates. In 3D, vortices 
are string-like objects, either forming closed loops (e.g., ring 
solitons) or terminating at two points on the condensate sur- 
face ID, |30[]. In this case, the interaction dynamics are in- 
stead characterized by string intercommutation or formation 
of smaller loops. 

Most of the BEC vortex literature has focused on vortices 
(or antivortices) in rotating traps, in which antivortices (vor- 
tices) are expelled from the condensate. The complex interac- 
tions described here instead rely on the presence of a popula- 
tion of vortices and antivortices and hence are best observed 
in a nonrotating system. Interfering condensates comprised of 
two, three, or more pieces will all generate the required con- 
ditions, provided they are nonrotating. Another example of a 
nonrotating system able to generate an equal population is the 
ring dark soliton system 02911 . Experimentally, phase contrast 
imaging with detuned light may leave the BEC sufficiently 
undisturbed to allow observation of these dynamics 113 ill . 

A vortex generates a local circulating velocity field and as- 
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FIG. 7: (Color online) Vortex interaction dynamics for the 120°, 
Io = 0.8 arrangement. In the lower frames vortices (dark) and an- 
tivortices (light) are shown, (a) Vortex trails corresponding to the 
range f=323-429 ms. (b) A propagating vortex-antivortex dipole 
(VD). f =430-471 ms. (c) Vortex-vortex-antivortex tripole rotat- 
ing through « 90°. /=47 1—517 ms. (d) Two VDs meet to form a 
quadrupole. VA partners are exchanged and the new VDs move off 
orthogonally to the original directions. /=526-537 ms. (e) A VD 
meets a lone vortex. The antivortex exchanges its vortex partner and 
the new VD moves off on a new trajectory. /=285-328 ms. (f) VA 
annihilation resulting from passing by a lone vortex. ?=342-366 ms. 

sociated phase gradient which falls off rapidly with distance. 
Another vortex or antivortex in this field experiences a force 
in the direction of flow 111 311 . Similarly, VDs travel at a ve- 
locity determined by their distance apart. In Fig. |7|a), the 
paths of several VDs, are shown over a 106 ms period. The 
dipoles in which the partners are more widely spaced have 
shorter paths, illustrating the predicted behavior (see also ac- 
companying movies [26]). In contrast, rotating vortex-vortex 
(or antivortex-antivortex) molecules are rarely and only fleet- 
ingly seen. Two equal-charge vortices can circulate around a 
common point, midway between them, like a facing pair of 
figure skaters. However, these structures have been shown to 
be unstable in radiative media 01 l32ll and would also be dis- 
turbed by the more mobile VDs, or prevented from rotating 
by (countering) field gradients due to the presence of nearby 
antivortices (vortices). 

In Fig. |7Jb), a VD approaches a lone, near-stationary an- 
tivortex, with which it couples to form the tripole shown in 
Fig. He). The net angular momentum of the tripole subse- 
quently causes it to rotate through approximately 90° before it 
interacts again H 1 5L 1 1 7f1 - In Fig.|3d), a rare event is shown, in 
which two counter-propagating VDs approach and momen- 
tarily meet, forming a quadrupole, before exchanging part- 
ners and moving off orthogonally to the original directions. 
In Fig. He), a VD approaches a stationary vortex, which is 
swapped for the partnered vortex. The new VD continues on 
a new trajectory. This case may be contrasted with the forma- 
tion of the tripole. In both cases, a VD approached a stationary 
(anti)vortex. However, the different interaction geometries re- 
sulted in the two different examples of dynamics shown. 

The spacing of VDs changes in response to local field gra- 
dient perturbations, causing them to slow, speed up, sepa- 
rate completely, or annihilate. An example of annihilation, 
promoted by the proximity of another vortex, is shown in 
Fig.|7tf). Following the annihilation, scattered remnant waves 
travel ahead of the event location, dissipating the residual VD 
kinetic energy. Such interactions of the vortex with field per- 
turbations and radiation of energy as waves are examples of 



vortex-sound interactions [1]. As already mentioned, vortex- 
vortex molecules are unstable, emitting spiral waves as they 
travel apart. Indeed, the sound waves emitted by all dynami- 
cal vortex interactions travel outwards and are reflected by the 
trap, recombining to form a chaotic fluctuating background 
that affects the ongoing condensate evolution. 

Parker et al. If33ll analyzed a BEC containing a lone vortex 
within a harmonic trap that was carefully modified to control 
the emission of spiral waves from a confinement region near 
the trap center. When spiral wave radiation was permitted to 
escape confinement and prevented from reinteracting with the 
vortex, the local energy was thus reduced, manifesting as a 
migration of the vortex away from the condensate center and 
its precession around the trap. Although our system is com- 
plex in comparison, we nevertheless observe sound emission 
from isolated vortices in the form of spiral waves. As our sys- 
tem evolves, a reduction in the number and associated energy 
of the vortices also occurs, consistent with the conversion of 
this energy to sound. 

Vortices residing in the outer parts of the trap spiral heli- 
cally about the center in a right-hand screw sense. Antivor- 
tices spiral in the opposite sense and are therefore likely to 
meet the aforementioned vortices. These often form VDs, 
which then move inward toward the trap center. 



VI. CONCLUSION 

A nonrotating pancake-shaped Bose-Einstein condensate 
(BEC) fragmented into three pieces, with a repulsive non- 
linearity, forms significant numbers of vortices and antivor- 
tices after merging. We have demonstrated, with numerical 
2D simulations of this system, that the vortex creation mecha- 
nism may be explained in terms of a linear theory of the inter- 
ference of expanding wave packets. This was contrasted with 
the vortex creation mechanism from a two-fragment BEC, in 
which dark stripe solitons decay into necklaces of vortices and 
antivortices; an intrinsically nonlinear physical process. With 
the three pieces separated sufficiently and arranged symmetri- 
cally, the formation of a distorted honeycomb lattice contain- 
ing equal numbers of vortices and antivortices was demon- 
strated and explained via the linear theory. Moreover, this 
theory shows that phase differences between the initially sep- 
arated pieces manifest as a global lattice translation. If al- 
lowed to expand in the absence of a trapping potential, the 
honeycomb lattice maintains its form as it expands. If instead 
the BEC evolves within a trap, the lattice was shown to melt, 
exhibiting a diversity of vortex interactions, including the for- 
mation of vortex-antivortex dipoles and other vortex clusters. 
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APPENDIX A: VORTEX LATTICE PRODUCTION 

In this Appendix we provide an analytical description of 
the vortices generated through linear superposition of three 
wavefunction fragments expanding from concentrated Gaus- 
sian distributions. The exposition here follows a similar ap- 
proach to that reported in 11211 which studied the vortices in the 
interference pattern from the Young's three-pinhole interfer- 
ometer. In that work, an analytical description of the far-field 
vortex locations was derived as a function of source arrange- 
ment for three equal-amplitude complex scalar waves repre- 
sented as either spherical waves or pinhole secondary sources. 
A representation in terms of a discrete parameter space arose, 
allowing estimates of the number of vortices — shown to relate 
to the information capacity of a beam — and the description of 
a natural coordinate system in terms of a family of hyperbolas. 

Consider three (2+l)D Gaussian wave packets, each of unit 
mass, whose position expectation values are stationary and 
centered about coordinates ri = 0, T2, r^. This is shown in the 
polar coordinate system schematic, Fig.[8ja). The normalized 
probability amplitude for a Gaussian centered at the origin, as 
a function of position r and time t is given by [34] 

w A n'^Ap/ti ( -(Ap/nf\r\ 2 \ 

*P(r, t) = t exp t , 

1 + i(Ap) 2 t/(mh) F \2[1 + i(Ap) 2 t/(mh)]J 

(Al) 

where m is the mass of an atom of the atomic species, Ap is 
the momentum uncertainty that defines the wave packet width, 
and r, = (r, f) is a (2+l)D position vector covering the xy- 
plane, whose origin lies at (r, t) = (0, 0). The total probability 
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amplitude arising from the three superposed Gaussians is then jj = a (t\ — 2r 3 rcos(8 - 63)) 



An 

b 3 = — + Inn, (A6b) 



T(r, t) 



Y n- ll2 Ap/n 
j-J 1 + i(Ap) 2 t/(mh) 



x exp 



f -(Ap/ft) 2 |r-r/ 
2[1 + i(Ap) 2 t/(mX)] 



+ 10 



(A2) 



where |r — r ,■] = yj(x- xj) 2 + (y - yj) 2 is the distance from the 
center of the j'th Gaussian to a given observation point (r, f) 
and <pj is the relative phase of the jth Gaussian. Separating 
the amplitude and phase terms of the wave field components 
Tj(r, t) = Aj(r, t) expfo/r, f)], we get 



Aj(r,t) 



and 



l/2 Ap/h 



1 + i(Ap) 2 t/(mn) 



exp 



/-(Ap) ; 



2 m 2 \r ■ 



\2[{Ap) 



> 4 t 2 + m 2 h 2 ] 



(Ap) 4 mt\r ■ 



2[(Ap) 4 t 2 n + m 2 h 3 ] 



(A3) 



(A4) 



Vortices lie at points of the phase ^(r, f), at which a line 
integral of V ± x along a closed path F about such a point eval- 
uates to a nonzero value. More formally § T dx = 2nn for 
some integer n + 0. At any point coinciding with a vortex 
core, the phasors pi,p2,P3 corresponding to the three wave 
packets, must sum to zero. Equivalently, Eq. (lA2b is equated 
to zero. With reference to Fig.[8jb), if the phasors are of equal 
length — corresponding to equal amplitude contributions from 
the three Gaussians — an equilateral triangle is formed, allow- 
ing the angles at the vertices to be specified simply. 

As the amplitude contributions from the Gaussians are not 
unconditionally equal we must make an appropriate approxi- 
mation to the amplitude term. By restricting consideration to 
some finite region of the xy-plane defined by |r - ry| < |r max - 
x$\, the exponential term in Eq. ( IA3b will be approximately 



unity provided its argument is small, or {Ap) t + m h » 
(Ap) 2 m 2 |r max - r/| 2 /2. For a given r max , this is always true af- 
ter sufficient time has elapsed. The amplitude term may now 
be factored out of the probability amplitude expression. Ex- 
panding the term |r-r ; | 2 = r 2 -2r-r ; +r i -r ; and remembering 
that ri = 0, the probability amplitude ^(r, t) from Eq. 
vanishes when 



1 + exp {/ |a [r\ - 2r 2 r cos dj + ^2]} 
+ exp [i [a (rj - 2r 3 r cos(6» -63)) + <fo]} = 0, 



(A5) 



where a = mht/[2(nt) 2 + 2m 2 (h/ Ap) 4 ]. The three summands 
are associated with the three phasors in Fig.[8jb). Whereas the 
first summand, 1, is uniquely identified with the horizontal 
phasor pi, the association of the sources with the other two 
phasors allows two permutations, corresponding to the two 
phasor diagrams; one for each vortex, vi and V2. In fact one 
will be a vortex and the other an antivortex. 

The arguments of the two exponentials in Eq. ( IA5t are 
denoted by y and rj, respectively. These phase angles are 
uniquely defined to within an integer multiple of 2n, so that: 



r 



; (r\ - 2r^r cos 6j + 1 



2 * „ 
h 2mn, 

3 



(A6a) 



where m and n are integers. The association of the m and 
n indices with the vertices is arbitrary. The choice is made 
here to match m with ri and n with r 3 . The alternative phasor 
association is thus 



£ = a [r\ - 2rir cos 0j + 1 



An „ 

h 2mn, 

3 



K = a{rl - 2r 3 rcos(6>- 6» 3 )) + (f> 3 = — + 2nn. (A7b) 



(A7a) 



Forming the fraction y/77 or (/k yields 



— ( cos 63 + tan tfsin #3 1 = — . 

r 2 \ ) rf-pM(m) 



r\-pN{n) 



(A8) 



where (3 = l/(3o) = 2[{ht) 2 + m 2 (h/Ap) 4 ]/(3mht). For the 
fraction y/rj, we get 



M(m) = 2n [1 + 3 (m - 4> 2 /2n)] , 
N(n) = 2n[2 + 3(n- fa/2n)] , 

and for the fraction £/k, we instead have 

M(m) = 2n[2 + 3(m- (p 2 /2n)] , 
AT(n) = 2^-[l +3(n-03/2n)]. 



(A9a) 



(A9b) 



Examination of these expressions reveals that a 2n change in 
relative phases (f>2 or (pj may be absorbed as an integer change 
in the associated parameter-space coordinate m or «, respec- 
tively. By allowing m and n to correspond to discrete real 
values instead of integer values, they may then also absorb 
fractional parts of the relative phase. Consequently, the vortex 
lattice may be continuously translated by a single lattice cell 
for each 2n change in the source phase. The position of any 
vortex in a cell coinciding with the BEC trap center is thus 
understood as resulting from a translation of the entire lattice. 
Isolating 6 in Eq. ( IA8b yields 



8 = arctan 



1 



r 3 -/3N(n)/r 3 
sin ft \r 2 -/3M(m)/r 2 



- COS O3 



(AlO) 



Finally, the expression for the radial coordinate r of the 
(m, n)th vortex core is obtained from Eq. ( IA6al ) or Eq. ( IA7al ) 
by applying the identity cos 2 # = 1/(1 + tan 2 #) and making 
use of Eq. dAlOb : 



1 



sin 2 6> 



pN{ri) 



[3M(m) 



COS #3 



+ 1. 



n \ r 2 

(All) 

Together with Eq. ( lAlOb . we have the vortex coordinates 
(r, 6, t), with the sign of the vortex charge — i.e. whether they 
describe the positions of vortices or antivortices — depending 
on the choice of Eqs. ( |A9at or dA9bl > and the value of 83. If 
83 e (0, 7r), vortices are indicated by Eq. dA9al > and antivor- 
tices by Eq. JA9bt . If 63 e (n,2n), the association is re- 
versed, with vortices indicated by Eq. (|A9b| ) and antivortices 
by Eq. dA9at . 
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Note that, in contrast with the case described in |[l2tl . in 
which the allowed range of integers (to, n) was restricted, 
here there are no restrictions and an infinite, uniform vortex- 
antivortex lattice is generated. However, the amplitude term 
Eq. (IA3b applies a Gaussian envelope to the probability den- 
sity I*?! 2 = |2j-=iA/'( r )| , effectively limiting the lattice. For 



sources arranged at the three corners of an equilateral trian- 
gle, the lattice has a symmetric honeycomb symmetry, with 
regular hexagonal cells. Changing the angle #3, or the side 
lengths r2 or r-}, distorts the cells and the lattice. 



